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ABSTRACT 

A numerical method Is developed to analyze the Invlscld flowfleld of a 
high speed Inlet by the solution of the Euler equations. The LU Implicit 
scheme In conjunction with adaptive dissipation proves to be an efficient and 
robust nonoscl 1 latory shock capturing technique for high Mach number flows as 
well as for transonic flows. 

INTRODUCTION 

Recent Interest in the aerospace plane and other hypersonic vehicles 
revitalized the research on high speed propulsion systems as well as hypersonic 
aerodynamics. In the design of supersonic and hypersonic propulsion systems 
the analysis of high speed flow past an Inlet plays a critical role. While 
there are half a dozen propulsion study concepts for high speed flight, many 
candidate concepts share a common Idea of combination of turboramjet engines 
for sub and supersonic flights and scramjet (supersonic combustion ramjet) 
engines for hypersonic flight. Resulting speeds of flows through the engines 
range widely from subsonic to hypersonic regimes. 

The analysis of hypersonic flows would require the full Navler-Stokes 
equations with slip effects and chemical reaction. It Is also Important to 
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understand the complex structure of shock waves. The turboramjet Inlet 
flowfleld Includes the Incoming supersonic flow deflected by oblique shock 
waves and the subsonic diffuser flow after the terminal normal shock wave 
while the scramjet Inlet flow Is characterized by strong oblique shock waves. 
The Euler equations which represent hyperbolic conservation law can be a useful 
testbed for developing and evaluating a shock capturing numerical algorithm. 

It has been a difficult task for computational aerodynamiclsts to capture 
nonosclllatory shock waves as a converged solution. Unbounded growth of 
spurious oscillations often resulted In numerical Instability. It is well 

known that upwind difference schemes can eliminate oscillations in the 

* 

neighborhood of shock waves at the expense of a substantial Increase of 

computational work. In parallel with the developments In upwind schemes it 

has been found that steady aerodynamic flows containing moderately strong shock 

waves can be quite well predicted by a central difference scheme augmented by 

a carefully controlled blend of first and third order dissipative terms. ^ 

In this paper the performance of adaptive dissipation Is demonstrated for 

strong oblique shock waves In high Mach number flows on a near-uniform mesh. 

Although a space marching method has been useful. It Is not well-posed 

when there Is upstream Influence through subsonic portions of the flowfleld 

such as a boundary layer. It also cannot handle the terminal shock and the 

subsonic diffuser flowfleld as well as the flow with streamwlse separation. 

Early time-integration codes for calculating the supersonic flow through an 
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Inlet used popular MacCormack schemes. ’ A disadvantage of these schemes 
for steady state calculation Is that the computed steady state depends on the 
time step. Another drawback of MacCormack's Implicit scheme Is the difficulty 
In treating boundary conditions. 

During the last decade, the Navler-Stokes equations have been the subject 
of exploratory Investigations aimed at establishing the feasibility of their 
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solution, but the methods so far developed have been too expensive to permit 

their use In a routine production mode. However, recent developments of modern 
4-6 

Implicit schemes In conjunction with multigrid methods are encouraging. The 
5 

authors developed an optimal dissipation model for the alternating direction 
Implicit (ADI) scheme and proved that the Improved ADI scheme Is Ideal for 
multigrid In two dimensions. Unfortunately, the ADI scheme In delta form to 
ensure the time step Independent solution has stability and convergence 
problems In three dimensions. Two new Implicit schemes which are 
unconditionally stable in any number of space dimensions were successfully 
developed by the authors^ recently. They are lower-upper (IU) Implicit 
scheme and LU-SSOR (symmetric successive over-relaxation) scheme. The LU 
Implicit scheme has been tried on an H-mesh In this work for high speed flow 
calculations; 

GOVERNING EQUATIONS 

The Euler equations are obtained from the Navler-Stokes equations by 
neglecting viscous terms. Let p, u, v, E, H, and p be the density, 

Cartesian velocity components, total energy, total enthalpy, and pressure, and 
let x and y be Cartesian coordinates. Then for a two-dimensional flow 
these equations can be written as 


aw aF aG n 

it *ii*i7'° 


(i) 


where W Is the vector of dependent variables, and F and G are convective 
flux vectors 

W = (p,pu,pv,pE) T 


2 T 

F = ( pu , pu + p,pvu,u(pE + p)) (2) 

2 T 

G = (pv ,puv ,pv + p,v(pE + p)) 

The pressure Is obtained from the equation of state 
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(3) 


P = p(y - l)|E - 1 (u 2 + v 2 )j . . 

These equations are to be solved for a steady state 3W/at = 0 where t 
denotes time. • •-.•••••••- 

' ' ' SEMI-DISCRETE FINITE VOLUME METHOO 

A convenient way to assure a steady state solution Independent of the time 
step Is to separate the space and time discretization procedures. In 
seml-di scre.te finite volume method one begins by applying a seml-dl scretl zation 
In which only the spatial derivatives are approximated.^ The use of a finite 
volume method for space discretization allows one to handle arbitrary 1 - 
geometries and helps one. to avoid problems with metric singularities tha.t ; are 
usually associated with finite difference methods. The scheme reduces to a. 
central difference scheme on a Cartesian grid, and is second order accurate in 
space provided that the mesh Is smooth enough. It also has the property that 
uniform flow Is an exact solution of the difference equations. 

; . NONLINEAR ADAPTIVE DISSIPATION 

In typical calculation of flow with discontinuities by a central 
difference scheme, wiggles appear In the neighborhood of shock waves where 
pressure gradient Is severe. In order. to suppress the tendency for spurious 
odd and even point oscillations, and to prevent unsightly overshoots near 
shock waves, the scheme Is augmented by artificial dissipative terms. The 
dissipative term, which Is constructed so that It Is of third order In smooth 
regions of the flow, Is explicitly added to the residual. For the density 
equation, for example, the dissipation has the form 

d 1 +1 /2 , j ' d 1-l/2,j * d 1,j+l/2 ' d 1,j-l/2 

where 


1+1/2, j 


- c (2) ( a 
~ c 1+l/2, j' p 1+l , j 


- p 


U 


j - 


(4) ' / o , o 

e 1+l/2,j (p 1+2,j ' Jp 1+l,j 3p 1,J 


p 1-l 

(4) 
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Let S be the cell area which Is equivalent to the Inverse of the determinant 
of transformation Jacobian. Both coefficients Include a normalizing factor 
S 1+l/2 j/At ;P ro P°r tJlonal t0 the length of the cell side, and j 

also made proportional to the normalized second difference of the pressure 


U 


huA : 2 p i- i : 


P H1J + 2P 1J + P 1-1J 


(5) 


In the adjacent cells. The third order terms provide background damping of 
high frequency modes. The first order terms are needed to control oscillations 
In the neighborhood of shock waves, and are turned on by sensing strong 
pressure gradients In the flow. The dissipative terms for the other equations 
are constructed from similar formulas with the exception of the energy equation 
where ; the differences are of pH rather than pE. The purpose of this Is to 
allow a steady state solution for which H remains constant. Increasing the . 
amount of artificial viscosity Improves the rate of convergence although too 
much dissipation can hurt It. However, It Is desirable to make the amount be 
as small as possible In order not to degrade the accuracy of solution. Typical 
amount of the third order terms Is almost negligible when compared to the 
physical viscosity. 

LU IMPLICIT SCHEME 
Let the Jacobian matrices be 

3G 


A - *£■ 
rt 3W 


B = 


3W 


( 6 ) 


and let the correction be 


6W = W 


n+1 


W 


here n denotes the time level 


The linearized Implicit scheme for a system of nonlinear hyperbolic 
equations such as the Euler equations can be formulated as^ 

(i + B fit (D A + D B)}iW + At R = 0 


(7) 
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where R Is the residual 

: R = D x F(W n ) D y G(w") - - ; , - , - 

Here D and D are central difference operators that approximate a/ax 
x y 

and 3/ay. ’ ' ' ' ' ' ‘ ' 

If B = 1/2 the scheme, remains second order accurate In time, for other 

values of 6, the time accuracy drops to first order. The unfactored Implicit 

scheme (Eq. (7)) produces a large block banded matrix which Is very costly : to 

Invert and requires hugh storage. One can solve the system Indirectly using a 

relaxation algorithm. Then it is desirable that the matrix should be 

diagonally dominant to meet a convergence criterion of a relaxation method. 

This can be achieved by flux splitting at the expense of a substantial Increase 

in the computational work. Moreover, it seems that second order flux splitting 

methods Tri conjunction with relaxation algorithms are either conditionally 

stable or slow. • 1 

The operation count cain be reduced by factorizing the operator of Eq. (7) 

approximately In various ways. The first way Is known as the ADI scheme.^ 

(I + B At 0 A)(I + B At D B) SW + At R = 0 (8) 

x y 

5 

Although the Introduction of optimal artificial dissipation makes the 
scheme be very desirable In two dimensions, the scheme In delta form is only 
conditionally stable In three dimensions. The ADI scheme Introduces error 

.V - - 3 

terms of order (At) In three dimensions which reduce the convergence rate. 

If one concerns about memory requirement, each factor can be split Into two 

O 

subfactors. If B = 1, the scheme becomes 

(I + At D~A + ) ( I + At D x A")(I + At D~B + ) ( I + At D^B'JAW + At R =0 (9) 

where D~ and D - are backward difference operators and D^ and D^ 
are forward difference operators. Each factor can be constructed using the 

9 

diagonally dominant ADI factorization. This scheme has six factors In three 
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dimensions and Introduces error terms of order (At)^ which reduce the 
convergence rate further.- ; 

While the ADI scheme has been valuable In two dimensions. Its Inherent 

limitations In three dimensions suggest an alternative approach. An 

unconditionally stable Implicit scheme which has error terms at most of order 
2 

(At) In any number of space dimensions can be derived by the LU 
factorization.® • 

(i + (3 At (D~A + + D"B + )} {i + 6 At (dV + dV))aW + At R = 0 (10) 

x y x y * 


Here, A + , A , B + , and B are constructed so that the eigenvalues of 
matrices are nonnegative and those of matrices are nonpositive. 


A+ = 2 (A * r A I} * A - = l (A - r A I) 

B + = \ (B * r B I), B' = \ (B - r B I) 


( 11 ) 


where 

r A > max( |X A I ) , r B >max(|x B |) (12) 

Here, X. and X D represent eigenvalues of Jacobian matrices. Equation (10) 

A o 

can be Inverted In two steps. The LU Implicit scheme needs the Inversion of 
sparse triangular matrices which can be done efficiently without using large 
storage. This scheme has only two factors In three dimensions. Other forms 
of factorization In conjunction with flux splitting can be found In Ref. 10. 

For example. 


(i + 6 At (D x A + + D y B)}’(l + (3 At DV AW 4- At R = 0) (13) 

This scheme requires a block-tr.ldlagonal Inversion In one direction. If one 
wants to Include thin-layer viscous terms In the Implicit operator, this scheme 
may be useful. However, It does not seem to be necessary to Insert viscous 
terms Into the operator when only the steady-state solution Is desired. 11 
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------- -- ---- - - RESULTS- --- ' 

In order to test the performance of the LU Implicit scheme for high speed 
flow calculations, a two-dimensional model problem was selected! Figure 1 
shows a typical hypersonic Inlet and Fig. 2 shows a 54 by 32 H-mesh for a 
schematic high speed Inlet. This mesh was used In all cases except the 
terminal shock wave problem where a 104 by 32 mesh was used for better 
resolution of the normal shock wave. However, all figures of convergence 
history are the results on the 54 by 32 mesh for comparison. The ramp angle 
is 9° and the shoulder angle Is set to 0.5° for the terminal shock wave 
problem. 

At the inflow boundary all the flow quantities are specified, and they 
are extrapolated from the Interior at the outflow boundary for supersonic 
outflow... For the terminal shock wave problem the pressure was prescribed at 
the outflow boundary. 

Four plots including Mach number contours, Mach number and pressure along 
the centerline, and the convergence history are shown In each set of Figs. 3 
to 6. Terminal shock wave problem with freestream Mach number 2 Is shown In 
Fig. 3. Figures 4 to 6 are for supersonic throughflows with freestream Mach 
numbers 5, 10, and 20, respectively. The pressure plots show the values of 
the pressure normalized by freestream static pressure so that the strengths of 
shock waves can be compared. Two Indicators In the convergence histories are 
the maximum and the average density residual In logarithm scale. 

As the figures show the wave structures In high Mach number flows 
Including oblique shock waves, reflected shock waves, expansion fans, and the 
Interaction of shock waves with expansion fans are successfully captured. 

These results clearly demonstrate the capability of the present numerical 
method for high speed flows. Figures 5 and 6 show that the location of the 
shock waves Is hardly changed as the Mach number Increases from 10 to 20. 
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However, the pressure plots show that the strengths of shock waves are quite 
different. 

The convergence histories show that the residuals drop linearly and 
continuously. These prove the efficiency of the present numerical method. 
However, as the Mach number Increases the convergence rate Is slowed down. 

This problem will be Investigated In the future. Another difficulty 
encountered In high Mach number flows Is the sensitivity to the way of starting 
the solution procedure. Sudden Introduction of the boundary condition to the 
freestream uniform flow Is likely to cause numerical instability In high Mach 
number flows. Gradual Increase of the time step is found to be effective to 
fix this problem. 

CONCLUSION 

The LU Implicit scheme combined with the nonlinear adaptive dissipation 
Is successfully developed as a robust and efficient shock capturing method for 
high Mach number Inlet flows. It seems to be possible to Improve the 
resolution and the accuracy of shock waves by using a total variation 
diminishing scheme. Extension to the Navier-Stokes equations Is desirable for 
more- accurate simulation of the flow through Integrated high speed propulsion 
system. 
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Figure 1. -A typical hypersonic inlet. 



Figure 2. - 54x32 H-mesh for a schematic high speed inlet. Ramp angle 
9° and shoulder angle OP or 0.5°. 


(a) Mach number contours. 
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(d) Convergence history. 

Figure 3. - Mach 2 inlet with the terminal shock wave. 
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(a) Mach number contours. 





Figure 4. - Mach 5 inlet. 




(a) Mach number contours. 
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(d) Convergence history. 
Figure 5. - Mach 10 inlet. 
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